First-principles study of epitaxial strain in perovskites 
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Using an extension of a first-principles method developed by King-Smith and Vanderbilt 
[Phys. Rev. B 49, 5828 (f994)], we investigate the effects of in-plane epitaxial strain on the ground- 
state structure and polarization of eight perovskite oxides: BaTiOs, SrTiOa, CaTiOa, KNbOa, 
NaNbC>3, PbTi03, PbZr03, and BaZrC>3. In addition, we investigate the effects of a nonzero 
normal stress. The results are shown to be useful in predicting the structure and polarization of 
perovskite oxide thin films and superlattices. 

PACS numbers: PACS: 77.55. +f, 77.80.Bh, 77.84.Dy, 81.05.Zx 



I. INTRODUCTION 

Ferroelectrics are insulating solids of technological im- 
portance because of their ability to maintain an electric 
polarization that can be reoriented by the application 
of an electric fieldpi This property lends itself to techno- 
logical applications including microelectronic devices and 
computer memories. Among ferroelectrics, perovskites 
constitute a subclass that has been of theoretical and 
experimental interest since the discovery in 1945 of its 
first member, barium titanate (BaTiOa). This interest is 
motivated in part by the relative simplicity of their cu- 
bic crystalline phase. For a perovskite of general formula 
ABO3, this structure contains cations A at the cube cor- 
ners, a cation B at the center of the cube, and oxygen 
atoms at the center of the cube faces forming a regular oc- 
tahedron, as depicted in Fig.^a). Typically, perovskites 
are found in this cubic paraelectric phase at high temper- 
ature; as the temperature is reduced, symmetry-lowering 
distortions to other phases, including ferroelectric ones, 
may occur. 

The electronics industry's demands for smaller com- 
ponents have made thin ferroelectric films the subject of 
recent attention^ 3 . Experimentally, it is found that the 




FIG. 1: (a) Ideal cubic perovskite structure for an ABO3 com- 
pound; the 2-polarized soft-mode atomic displacements are 
indicated by arrows, (b) Epitaxial paraelectric (or p) phase, 
in which the atoms are constrained in the xy plane due to the 
presence of the substrate. 



properties of ferroelectrics in thin-film form generally dif- 
fer significantly from those in the bulk. While many fac- 
tors are expected to contribute to these differences, it has 
been shown that the properties of perovskite thin films 
are strongly influenced by the magnitude of the epitax- 
ial strain resulting from lattice-matching the film to the 
substrate, known as misfit strain or epitaxial strain. 

Previous theoretical studies have isolated the effects of 
epitaxial strain on the structure and properties of films by 
imposing the epitaxial constraint on the in-plane lattice 
vectors of a periodic bulk sample. Using a phenomeno- 
logical Landau-Devonshire model, Pertsev, Zembilgotov 
and Tagantsev^ introduced the concept of mapping the 
equilibrium structure of a ferroelectric perovskite mate- 
rial versus temperature and misfit strain, thus produc- 
ing a phase diagram of the observable epitaxial phases. 
Given the importance of strain in determining the prop- 
erties of these films, these diagrams have proven to be 
of enormous interest to experimentalists seeking to in- 
terpret the results of experiments on epitaxial thin films 
and heterostructures. This phcnomcnological approach 
should give excellent results in the temperature/strain 
regime in which the model parameters were fitted (usu- 
ally near the bulk ferroelectric transition) but will gener- 
ally be less accurate when extrapolated to other regimes. 
In particular, Figure 1 of Ref. |5| shows that two differ- 
ent sets of parameters can give two quite different phase 
diagrams. Furthermore, it is only possible to study ma- 
terials for which all the needed experimental information 
is available. 

In previous workj^Sii we have examined the effects 
of epitaxial strain with an analogous, but fully first- 
principles, approach. Specifically, we presented density- 
functional theory (DFT) 8 calculations for the structure 
and properties of BaTi0 3 , PbTi0 3 and SrTi0 3 with 
varying in-plane strain, fully relaxing all structural de- 
grees of freedom consistent with uniform distortions (that 
is, retaining the five-atom unit cell). From this, we 
obtained zero-temperature phase diagrams that com- 
plement the phcnomcnological results of Pertsev and 
coworkers. In this paper, we show that these phase 
diagrams can be reproduced using the first-principles 
energy parameterization of King-Smith and Vanderbilt, 
and give results for an additional five perovskite oxides: 
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CaTi0 3 , KNb0 3 , NaNb0 3 , PbZr0 3 , and BaZrOg. This 
approach greatly decreases the computational effort in- 
volved in computing the phase diagram and readily al- 
lows the inclusion of nonzero external stressiS Moreover, 
the parameters can be fully specified in a compact ta- 
ble, and the functional form of the energy is suitable for 
analytical computations and conceptual interpretation, 
leading to a classification of possible stress-strain phase 
diagrams. These results can be used for predictions of 
the structure of epitaxially strained thin films grown on 
substrates with square symmetry, and for the design of 
novel perovskite strained-layer superlattices with two or 
more components. 

The paper is organized as follows. In Sec. [H] we de- 
scribe the extension of the KSV method to study the ef- 
fects of the epitaxial strain constraint and external stress, 
and the first-principles calculations used to obtain the 
KSV parameters. Section IIIII presents the results for 
the eight perovskite oxides considered, including the se- 
quence of phase transitions with varying misfit strain at 
zero stress, and phase diagrams in which we show the 
most stable phase for given misfit strain and external 
stress. In Sec. IIVI we review the approximations and 
discuss how to use the present results for the prediction 
of the structures and polarization of thin films and su- 
perlattices. Finally, in Sec.^we present our conclusions. 



II. METHOD 
A. Formalism 

The starting point of this analysis is the parameter- 
ized total-energy expression presented by King-Smith 
and Vanderbilt in Ref. |3 This is a Taylor expansion 
around the cubic perovskite structure in terms of the six 
independent components r)i of the strain tensor (i is a 
Voigt index, i =1-6) and the three Cartesian components 
u a (a = x, y, z) describing the amplitude of the soft mode 
defined by the pattern of eigen-displacements associated 
with the smallest eigenvalue of the (zone-center) force- 
constant matrix. The arrows in Fig. ^ indicate a typical 
displacement pattern associated with this mode, greatly 
magnified to allow its visualization. 

The contributions to the energy (per unit cell) can di- 
vided into terms arising from pure strain and from pure 
soft-mode amplitude, and an interaction term, 

E = E^({th}) + E s ° a ({u a }) + E^({ Vl }, {u a }), (1) 

with the zero of the energy corresponding to the cubic 
structure. For crystals with cubic symmetry the strain 
energy is given, correct to second order in the strains, by 



E^({ m }) 



1 



+B 12 {mm + mm + mm) 

+ ^B 44 ( v 2 + t 1 2 + t 1 2 ), 



where Bn, B12, and -B44 are related to the elastic con- 
stants of the crystal by factors of the cell volume. The 
soft-mode energy given in Ref. ^3 contains terms up to 
fourth-order in the soft-mode amplitude, 

E solt {{u a }) = nu 2 + au 4 + i{u 2 x u 2 y + u 2 y u\ + u 2 z u 2 x ), (3) 

where u 2 — u 2 + u 2 + u 2 , k is twice the soft- mode eigen- 
value, and a and 7 are the two independent symmetry- 
allowed fourth-order coefficients describing the cubic 
anisotropy. Finally, the interaction between the strains 
and the soft-mode amplitude is given by 



E int ({ Vl },{u a }) = -B lxx ( Vl u 2 x + mu 2 y + Vsu 2 z ) 



+ 7, B iyy[m{u y + u z ) + T] 2 (u z + u x ) 

+m(ul + ul)] 
+B iyz {mu y u z + m u zu x + m u xu y ), (4) 

where B\ XXl B\ yy , and B 4yz are the phonon-strain in- 
teraction coefficients. All the coefficients in these three 
parts of the total-energy expression can be obtained from 
first-principles calculations on a series of distorted struc- 
tures as described in Ref. [10| and in the next subsection. 

In this paper we will be concerned with the effects of 
strain on a film grown epitaxially on a substrate with 
square symmetry. The epitaxial strain constraint im- 
posed by the substrate is 



m = m = v, 
m = o, 



(5) 
(G) 



where fj is the misfit strain between the minimum-energy 
cubic structure of the film material and the substrate. 

In the case of epitaxy, where strain elements 771, 772 
and 776 are constrained while the others are not, it is 
useful to introduce a mixed stress-strain elastic enthalpy 
G = E — 03773 — 174774 — 05775 whose natural variables 
are u x , u y , u z , 771 , 772 , ?76 , 03 > °4 an d 05 . Specializing to our 
case in which 771 = 772 = m V6 = 0, an d assuming that 
the shear stresses 04 and 05 vanish, we define an effective 
elastic enthalpy given by 



G = E - CT3773 



(7) 



whose natural variables are u x ,u yi u Zl fj and er 3 . Using 
Eqs. ill'l II) and minimizing Eq. J7J with respect to 773, 774 
and 775 yields 

% = -5— [03 - 273i277 

-\Bi xx u 2 z --B lyy {u 2 x + u 2 y )] (8) 



B_ 



(2) 
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(9) 
(10) 
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Substituting these expressions back into Eq. J7J), we ex- 
press G in terms of its natural variables as 



G 



(A vri r] 2 
+(Bfjfj- 



CTCT 



-Du 



XI) 



-Eui 



■ B)u 2 xy 
C)ul 



-Hu xy sin 2 9 cos 2 9. 



11) 



Here we have simplified the notation by replacing 03 by 
a. Also, the two soft-mode amplitude components are 
represented in polar coordinates as 



u „ 



, cos 9, 
, sin 9. 



(12) 
(13) 



The coefficients in G are expressed in terms of the KSV 
parameters as follows: 



B 2 

Afjfj = Bn + B12 — 2 — — , 
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(14) 

(15) 

(16) 

(17) 

(18) 
(19) 
(20) 

(21) 
(22) 

(23) 

(24) 

(25) 
(26) 



For a given set of coefficients in the potential G of 
Eq. 1)11(1. we can predict the phase diagram as a function 
of misfit strain fj and the normal external stress a by 
minimizing G to find the values of the ground-state soft- 
mode amplitude components. For a fourth order theory, 
like the present KSV expression, the entire optimization 
process can be done analytically, since it is possible to 
compute first and second derivatives of G and to do a 
stability analysis of the various possible phases, classified 
by the nature of the minimum-energy soft-mode vector. 
For example, a paraelectric p phase similar to the cubic 



TABLE I: Summary of epitaxial perovskite phases. Columns 
list, respectively, the phase, space group, and symmetry of 
the soft-mode amplitude components. 



Phase 



SG 



SMA components 



c 
a 

aa 
ac 
r 



PAmmm 
PAmm 
Pmm2 
Amm2 

Pm 

Cm 



U X = Uy = U Z — 
«I = Uy = 0, U Z 

U x ^ 0, Uy = u z = 
u x = u y / 0, u z — 

U x / 0, Uy = 0, U z 

U X = Uy / 0, u z 



phase can appear if the potential is minimized for u = 0, 
but relaxation along z will occur, making the cell tetrag- 
onal, as shown in Fig^b). The classification, following 
Pertsev and coworkers, 4 is given in table [I] Expressions 
for the elastic enthalpy of a given phase as a function 
of misfit strain can be obtained by minimizing Eq. 1)11(1 
with the appropriate constraint on u. For example, the 
elastic enthalpies for the p, c, a and aa phases are 



(27) 
(28) 

(29) 

(30) 
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The process of generating the stress-strain diagrams 
is thus extremely rapid in comparison with a full DFT 
analysis, once the KSV parameters have been obtained 
from first principles calculations as described in the next 
section. Justification of the approximations involved has 
been presented in Ref. 0, and will be discussed further 
in Sec.HVI 



B. First-principles calculations of the coefficients 

In Table V of their paper^ King-Smith and Van- 
derbilt report the computed coefficients to be used in 
their model. We have now repeated calculations anal- 
ogous to theirs, taking advantage of the increase in 
computational power that has taken place during the 
last ten years to push the boundaries of the numerical 
approximations that control the accuracy of the first- 
principles calculations. We used the same local-density 
approximationii*i£ DFT methodology and the same ul- 
trasoft pseudopotentials that they used>i£ In our case, 
the plane-wave kinetic-energy cutoff has been raised to 
50 Ry, and the Monkhorst-Pack 14 k-point mesh is finer, 
containing 8x8x8 points. 

The lattice parameters and the soft-mode eigenvector 
components calculated for the eight perovskites under 
study are reported in Table ^ These values are quite 
similar to those reported in Tables IV and VII of Ref. 0, 
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TABLE II: Lattice parameters (in bohr) and soft-mode eigen- 
vector components (normalized to unity) for eight perovskites, 
calculated using DFT. 





ao 


£a 


£b 


£Ol,2 


£o 3 


BaTi0 3 


7.448 


0.184 


0.774 


-0.218 


-0.522 


SrTi0 3 


7.285 


0.490 


0.596 


-0.408 


-0.270 


CaTi0 3 


7.201 


0.678 


0.363 


-0.435 


-0.171 


KNbOa 


7.470 


0.195 


0.796 


-0.325 


-0.341 


NaNbOa 


7.395 


0.427 


0.638 


-0.428 


-0.209 


PbTi0 3 


7.337 


0.597 


0.487 


-0.411 


-0.262 


PbZr0 3 


7.762 


0.769 


0.152 


-0.438 


-0.044 


BaZr0 3 


7.846 


0.703 


0.275 


-0.462 


-0.054 



where a comparison with experimental and other theo- 
retical data can also be found. 

Following the prescription given by King-Smith and 
Vanderbilt^S we found the updated KSV parameters dis- 
played in Table IIIII The two final columns of this ta- 
ble contain the values of a' and 7', the coupling con- 
stants whose values determine the symmetry of the low- 
temperature phase for bulk perovskites^ (see Ref. IT(i[) . 
Most of the parameters are only slightly different from 
those given in Table V of Ref. llOl and we have found 
that the predictions that both sets give for the epitax- 
ial films are qualitatively the same (note, however, that 
cubic SrTi0 3 changes from being marginally unstable to 
marginally stable) . The differences are mainly related to 
the use of finer grids in our case when performing the fast 
Fourier transforms required in the calculations. On the 
other hand, the improvements in plane-wave cutoff and 
k-point mesh have a quite small impact on the values of 
the coefficients. 

For our discussion of the effects of epitaxial strain on 
the polarization, it is useful to have an expansion not just 
of the energy, but also of the polarization, in terms of the 
soft-mode amplitude. The linearized expression 

P z = ^z*u z (31) 

is adequate for most purposes. Here e is the absolute 
value of the electron charge, f2 is the unit cell volume of 
the perovskite, and z* is the Born effective charge of the 
soft mode, given in terms of the soft-mode eigenvectors 
of Table [FT] by 

z* = ZXU + Z&b + 2Z* 0i J 0l . 2 + Z^o 3 - (32) 

We take the Born effective charges Z* for each atom 
to have their values in the cubic structure as calcu- 
lated by Zhong, King-Smith, and Vanderbilt. 16 Using 
these together with the eigenvectors reported in Table 
|n]we obtain the following values for z*\ 9.94 (BaTi0 3 ), 
8.65 (SrTi0 3 ), 7.03 (CaTi0 3 ), 11.06 (KNb0 3 ), 9.14 
(NaNb0 3 ), 9.40 (PbTi0 3 ), 6.29 (PbZr0 3 ), and 5.69 
(BaZr0 3 ). This approximate expression neglects the pos- 
sible dependence of the Born effective charges upon the 
strain. 



III. RESULTS 
A. Calculations at zero external stress 

We first consider the case in which the external perpen- 
dicular stress a vanishes. Table llVl shows the sequence of 
transitions that occurs for each of the eight perovskites 
studied as the misfit strain increases, and the values of 
strain at which the transitions occur. 

The observed sequences of phases can be understood 
with the help of Fig. |3 which illustrates the types of 
elastic enthalpy behaviors that we observe for the ma- 
terials considered. At the strains at which a transition 
occurs from one phase to another, the energy curves join 
smoothly, indicating that these transitions are of second 
order. It can be shown analytically that this is indeed 
the case for the KSV model, and that at the symmetry- 
breaking transitions (c-r, aa-r, a-ac, p-c, p-a, or p-aa) 
the higher-symmetry phase becomes unstable. For all 
compounds considered, we see that for sufficiently high 
compressive strains, the lowest energy phase is always the 
c phase, in which the atomic displacements and therefore 
the polarization point in the [001] direction, perpendicu- 
lar to the substrate. On the other hand, for sufficiently 
high tensile strains, we obtain a phase in which the po- 
larization lies in the substrate plane, pointing along [100] 
(a phase) for for BaZr0 3 or along [110] (aa phase) for the 
rest of the compounds. The in-plane orientation is deter- 
mined by the sign of H (i.e., of 7), which is positive for 
BaZr0 3 and negative for the seven other compounds. In 
the intermediate strain regime, three different behaviors 
are found. For BaTi0 3 , KNb0 3 , NaNb0 3 , and PbZr0 3 , 
an r phase appears between the c and aa phases, as in 
Fig. |2t a) . The fact that these perovskites crystallize in 
the r phase at low absolute values of strain is not sur- 
prising, since this phase is the most similar to the rhom- 
bohcdral phase that they adopt as the bulk ground state 
according to the KSV theory (see Table VI of Ref. 
For SrTi0 3 and BaZr0 3 it is the paraelectric p phase 
that appears at intermediate strains, as in Fig.|2Ib). For 
these two compounds, no r or ac phase appears. Instead, 
the polarization along [001] continuously goes to zero as 
the strain becomes less compressive, and only reappears 
in the xy plane once the tensile strain reaches some given 
value, continuously growing from then on. The epitax- 
ial paraelectric p phase is the analog of the bulk cubic 
phase (see Fig.^b)), which is the ground state predicted 
by the KSV theory with the parameters of Table IIIII for 
both bulk SrTi0 3 and bulk BaZr0 3 . Finally, for CaTi0 3 
and PbTi0 3 yet another behavior is obtained, as shown 
in Fig. |5fc). At intermediate strains, the rhombohedral 
phase is the lowest energy single phase. However, partly 
because of its inverted-parabola elastic enthalpy curve, 
the common tangent line between the c and aa phases 
yields a lower energy in the intermediate strain regime, 
and thus a mixed phase of c and aa domains is expected. 

The strain-induced phase transitions found using this 
approach compare well with the full DFT results previ- 
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TABLE III: Coefficients of energy expansion, Eqs. Ill It . for eight perovskites, in atomic units. Coefficients a' and 7' are not 
directiy reievant for this work, but are inciuded for compieteness (see Ref. 





B11 


B12 


Baa 


B\xx 


Blyy 


Biy Z 


K 


a 


7 


a 


7 


rial 1U3 




Loll 


1. 1 y 


— Z.Zl 


ft 1(1 

— u.iy 


n 1 n 
— U.1U 


— U.UMo 


U.o4 / 


O A QQ 

— U.400 


n 1 no 


n 1 on 


SrTi0 3 


4.97 


1.50 


1.54 


-1.28 


0.06 


-0.12 


0.0010 


0.124 


-0.161 


0.074 


-0.036 


CaTiOa 


5.34 


1.24 


1.41 


-0.65 


0.09 


-0.11 


-0.0070 


0.031 


-0.017 


0.019 


0.012 


KNbOs 


6.66 


1.03 


1.46 


-2.91 


0.37 


-0.03 


-0.0155 


0.439 


-0.655 


0.257 


-0.178 


NaNbOs 


6.15 


1.18 


0.98 


-1.71 


0.51 


-0.02 


-0.0130 


0.209 


-0.298 


0.124 


-0.050 


PbTiOs 


4.58 


1.86 


1.40 


-0.79 


-0.05 


-0.07 


-0.0117 


0.052 


-0.020 


0.031 


0.029 


PbZrOa 


5.34 


1.64 


0.93 


-0.20 


0.06 


0.01 


-0.0185 


0.013 


-0.013 


0.011 


-0.008 


BaZrQ 3 


5.69 


1.46 


1.45 


-0.44 


0.05 


-0.12 


0.0087 


0.013 


0.003 


0.008 


0.012 



TABLE IV: The sequence of epitaxially-induced phase tran- 
sitions and the vaiues of strain fjx and 7711 at the boundary 
of the c phase and of the aa (or a) phase, respectiveiy. An 
asterisk denotes the strain regime where formation of mixed 
domains of c and aa phases couid be favorable. 





Sequence 


n±(w-' A ) 


7711 (10- 3 ) 


BaTiOs 


c-r-aa 


-5.89 


7.59 


SrTi0 3 


c-p-aa 


-2.24 


1.59 


CaTiOa 


c-*-aa 


-2.30 


5.35 


KNbOs 


c-r-aa 


-4.80 


5.49 


NaNbOs 


c-r-aa 


-5.52 


4.13 


PbTiOs 


c-*-aa 


-3.00 


8.42 


PbZrOs 


c-r-aa 


-52.06 


30.42 


BaZrOs 


c-p-a 


-53.41 


41.86 



ously reported for BaTi0 3 £ PbTi0 3 £ and SrTi0 3 i We 
look now in detail at this comparison for BaTi0 3 . Figure 
13 shows the energy curves of the various phases as pre- 
dicted by the KSV theory (right panel), to be compared 
with the full DFT results (left panel). The agreement 
between the two sets of results is very good, with the 
small differences present arising from two sources. First, 
the first-principles calculations in Ref. [5| were performed 
using the projector-augmented wave method^ while the 
first-principles calculations used to obtain the KSV coeffi- 
cients in the present work were performed using ultrasoft 
pscudopotentials, 13 Second, there are the intrinsic errors 
due to the use of a Taylor expansion described in the pre- 
vious section, which are expected to grow as the strain 
and soft mode magnitudes increase. 

Figure 0] shows the displacements of the atoms from 
their centrosymmetric perovskite positions as strain 
varies. Again, the agreement of the KSV results (right 
panel) with the full DFT results (left panel) is very good. 
In particular, the square root behavior predicted by the 
KSV theory is exhibited by the more exact DFT calcu- 
lations. As the in-plane strain increases, we observe a 
second-order phase transition (c-r), and while the mag- 
nitude of the atomic displacements continues to diminish 
along [001], the displacements in the xy plane begin to 
grow. With increasing tensile strain, the displacements 
along [001] vanish at the r-aa transition, while the dis- 




Misfit Strain Misfit Strain 



FIG. 2: Sketches of different behaviors found for the elastic 
enthalpy curves of the most stable phases. Solid circles show 
where two parabolas meet. Empty circles show the points 
where a tie line meets a parabola. The curves that meet at 
the circles do so with equal first derivatives. 

placements in the xy plane continue to grow smoothly. 
In this way, we see that the polarization vector contin- 
uously rotates in going from the c phase through the r 
phase to the aa phase. A quantitative limitation of us- 
ing a single misfit-strain-independent local mode in the 
KSV model is shown here in the form of an artificial con- 
straint equating the A y (Oi) and Aa;(0 3 ) displacements. 
This constraint is removed when full DFT calculations 
are performed and the atoms are free to relax within the 
given space group, but the magnitude of these displace- 
ments is not very different from that obtained in the KSV 
theory. 

Figure [S] shows the values of P z for the various com- 
pounds as a function of misfit strain. The square-root 
singularity at P z =0 corresponds to an r-aa or c-p transi- 
tion. The slope discontinuity visible at finite P z in some 
curves corresponds to the c-r transition of Fig.^a), while 
the termination of the curves at finite P z for CaTi0 3 
and PbTi0 3 corresponds to the encounter with the tie 
line in Fig. 01b). The general trend, of course, is a 
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Direct DFT Calculations 



KSV Model 




-15-10 



Misfit Strain (10 



FIG. 3: Comparison of BaTiOa energy curves for the 
six epitaxial phases studied, as obtained from direct DFT 
calculations^ (left), and from the KSV theory (right). En- 
ergies are relative to the paraelectric structure at zero misfit 
strain. The lines in the left panel and the symbols in the right 
panel are provided as guides to the eye. 




7.3 7.4 7.5 

In-plane lattice constant (bohr) 

FIG. 5: Value of the polarization along the z axis for different 
perovskites. The continuous thick part of each curve indicates 
easily achievable misfit strain conditions, with misfit strain 
values between —20 x 10~ 3 and 20 x 10~ s , while the thin 
dashed part of each curve corresponds to larger strains. The 
curves for CaTiOs and PbZr03 fall outside the plotted region. 



B. Stress-strain phase diagrams 
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FIG. 4: Comparison of BaTiOa atomic displacements for the 
most stable phase at each given value of strain, as obtained 
from direct DFT calculations- (left), and from our KSV model 
(right). A z (Ti) indicates the displacement of the Ti atom 
along the z direction, etc. Symmetry implies that A^(Ti) = 
A K (Ti) , A*(0 2 ) = A„(Oi), A„(0 2 ) = A^Oi), A z (O a ) = 
Az(Oi), and A y (0 3 ) = A x (0 3 ). The lines in the left panel 
are guides to the eye. 



strong increase of P z with compressive strain, with con- 
siderable enhancement possible over the zero misfit strain 
value. Interestingly, according to this plot CaTiOa would 
have a large P z if the zone boundary distortions that are 
present in its actual ground-state crystal structure were 
suppressed. This result is supported by recent full DFT 
calculations i& 



We now consider application of a nonzero normal stress 
a. The stress-strain phase diagrams obtained for each of 
the eight perovskites are shown in Fig. H3 

All eight diagrams show a universal topology with 
straight-line phase boundaries meeting at a single cross- 
ing point. The perfect linearity of the boundaries is an ar- 
tifact of the truncation of our energy expansion to fourth 
order, but the presence of a single crossing point is robust 
against the introduction of small higher-order terms. The 
crossing point strain fj x and stress tr x can be connected 
with the critical strain and stress in the isotropic cubic 
perovskite at which a structural instability first occurs 
as the pressure is reduced or made more negative. There 
are two main variants of the diagrams in Fig. [HI the first 
with four phases (c, r, p, and aa/a), and the second with 
three phases (c, p, aa/a and a mixed phase region). The 
varied behavior of the zero-stress diagrams of Fig. [21 can 
now be interpreted, in the context of Fig. [BJ as reflecting 
whether the zero-stress axis lies above, or below, a x ■ 

More specifically, the coordinates of the crossing point 
can be expressed as functions of the total-energy coeffi- 
cients as follows: 



?/x = 



B a C — BC a 

BffC' a — B a Cfj 

BCfj — Bf t C 

BfjCa ~ BfjCfj 



2 k 



fjx{Bu 



2B\yy 

-2Bi 2 ). 



(33) 
(34) 



The strain and stress values at the crossing point are 
rather modest, with the single exception of PbZrC>3. 
PbZrC>3 has the lowest soft-mode eigenvalue, and there- 
fore the most negative value of k (Table llllf) . Its x- 
polarized soft-mode eigenvalue also changes more slowly 
on application of an 771 strain, resulting in a less negative 
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FIG. 6: External stress versus misfit strain phase diagrams for 
eight epitaxially strained perovskites. Straight lines represent 
second-order phase transitions. Shadowed areas indicate the 
presence of c and aa (or a) domains. 



value of B\ xx . These values of n and B\ xx result 

and c7 x being about an order of magnitude more negative 

for PbZrC>3 than for the other seven compounds. 

For external stresses below cr x , the behavior is sim- 
ilar in all eight compounds, showing a c-p-aa (c-p-a 
for BaZrOs) sequence of second-order phase transitions, 
with the elastic enthalpies of the phases behaving as in 
Fig. 01b). In this fourth-order KSV theory, the phase 



C_ 
B 

~bZ 



B,-i 
B, 



(35) 
(36) 



Rewriting these expressions as functions of the funda- 
mental coefficients of Table IIIII we see that the value of 



B\ V y plays a central role. For materials like the ones we 
are studying, for which S u > 0, J5 12 > 0, and B lxx < 0, 
the slope of the c-p transition line is found to be positive 
if 

|H > **L. (37) 

For all eight perovskites this slope is indeed positive, 
since the left hand side of the inequality is positive, and 
Biyy is positive or only slightly negative but much smaller 
in magnitude than Bi xx (the latter being the case for 
BaTiOa and PbTiOs). The slope of the p-aa transition 
is positive if 



B12 I 
B u > 2 



1 



Bl X3 

~B~ 



(38) 



where the inequality is not satisfied for BaTiC>3 and 
PbTiC>3, which have negative values for B\ yy . There- 
fore, for these two perovskites a transition from the aa 
to the paraelectric p phase is induced by applying a suf- 
ficiently high external tensile stress at fixed misfit strain, 
while for the others the transition would be from the p 
to the aa (or a, in the case of BaZrOs) phase. 

For external stresses above tr x , two kinds of behav- 
iors are found. Five of the perovskites (BaTiC>3, SrTiC>3, 
KNbC>3, NaNbC>3, and PbZrOs) show a c-r-aa sequence 
of second-order phase transitions under these conditions, 
with the energies of the phases behaving as in Fig. OJa) . 
In this case, the phase boundaries are straight lines of 
the form 

2BE - CF 



2B„E - C a F 



CfjF 



2B a E — C a F 
2C(D + if/4) - BF 
2C a (D + H / 4) — B a F 
2C n (D + H/4)-B n F 



>1- 



(39) 



(40) 



2C a (D + H/A)-B a F' 

For all five compounds, the slopes of both boundaries are 
positive. 

The other behavior observed for external stress above 
<7x is one in which c and aa (or a, in the case of BaZrOs) 
domains are expected. The energy curves behave either 
as shown in Fig. |2Ic) (for CaTiC>3 and PbTiOs) or as 
shown in Fig.^d) (for BaZrOa). However, in the inter- 
mediate region, instead of a uniform phase, the system is 
expected to break into domains, as explained in the pre- 
vious section for PbTiC>3 and CaTiC>3 at zero external 
stress. 
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IV. DISCUSSION 

In this section, we first consider in detail the several 
approximations that are responsible for the ease and sim- 
plicity with which we can generate stress-strain phase 
diagrams for perovskites. We then discuss and give ex- 
amples of the applicability of the theory to realistic ex- 
perimental studies of perovskite films and superlattices. 

Within the KSV theory the thermodynamical poten- 
tial is expanded as a Taylor series in strain and soft-mode 
amplitude, where the reference used is the perfect cubic 
perovskite of Fig. Qfa). The truncation at low order in 
the variables of the expansion means that the expansion 
decreases in accuracy for large distortions. As relevant 
misfit strains are generally rather small (less than 2%), 
this appears not to have significant implications. In addi- 
tion, in the present form of the theory, the modes selected 
for the expansion allow only phases that involve five-atom 
unit cells to be considered. In particular, we do not take 
into account the possibility of cell-doubling oxygen octa- 
hedra rotations, which have been shown to be important 
in SrTi0 3 ,i^a CaTi0 3 ^ and PbZr0 3 ^i In principle, 
such rotations could condense in the other compounds 
under high enough misfit strains, but for BaTi0 3 it has 
been shown^ that this does not occur until one reaches 
experimentally irrelevant strains, and we expect that this 
will also be the case for most of these other compounds. 
Second, our calculations are done at zero temperature. 
Extending them to finite temperatures could in principle 
be done using an effective Hamiltonian method as we did 
in Ref. 5 for BaTi0 3 . However, this would involve design- 
ing for each perovskite an effective Hamiltonian along the 
lines described by Zhong, Vandcrbilt, and Rabe^ and is 
left for future work. Third, our model does not include 
the small effect of the zero-point motion of the ions (see 
Ref. |2i| for a discussion). Finally, our theory relies on 
the LDA to compute the exchange and correlation terms 
in DFT. This introduces small systematic errors in the 
calculation, the most important of which is probably the 
error in the equilibrium lattice constant. However, such 
errors are well understood and well characterized in per- 
ovskites, and tend to be similar for different materials of 
this class, so that there is a tendency for cancellation of 
errors in relative quantities such as misfit strains (see, for 
example, Ref. ITot). 

The phenomenological Landau-Devonshire approach 4 
also requires approximations, which we review here for 
the purposes of comparison. Its starting thermodynam- 
ical potential is the bulk free energy expanded in polar- 
ization and stress, with linear temperature dependence 
in selected coefficients. Sixth order terms are needed as 
their importance increases at finite temperature^ The 
reference used is the paraelectric cubic perovskite phase 
at the bulk critical temperature T c , and the parameters 
are fit to reproduce experimental observations of the be- 
havior near the bulk ferroelectric transition. For the epi- 
taxial strain dependence, a Legendre transformation is 
then made to obtain the potential as a function of po- 



larization and misfit strain. With parameters extrapo- 
lated to zero temperature, this can be compared to the 
potential in the present work using the linear relation 
between the polarization and the soft mode amplitude 
(|31|l . Due to the way in which the parameters are fit, the 
Landau-Devonshire potential will give its most accurate 
results for small misfit strains and temperatures near the 
bulk T c , while the first-principles potential will be more 
reliable for the zero-temperature misfit-strain phase dia- 
gram. 

We now turn to the applicability of our theory to the 
prediction and understanding of properties of experimen- 
tally relevant systems. We first discuss the case of thin 
films, and then consider the case of strained-layer super- 
lattices. 

The theory presented here can be used directly to pre- 
dict the structure and polarization of a single-domain 
perovskite-oxide thin film grown on a substrate with 
square-lattice symmetry. The effects of epitaxial strain 
will be most evident in films coherent with the substrate. 
In equilibrium, coherent epitaxial growth is possible up 
to a certain critical thickness, which depends upon the 
misfit between film and substrate materials as well as 
upon temperature and other growth conditions. Us- 
ing low-temperature synthesis, coherence can be main- 
tained far beyond the equilibrium critical thickness, as 
has been shown for BaTi0 3 grown on GdSc0 3 (—1.0% 
misfit strain)^ Under such coherent conditions, the full 
misfit should be used as the input for our theory. For ex- 
ample, in the case of BaTi0 3 on GdSc0 3 , Fig.^shows a 
predicted enhancement of P z from 0.21 to 0.31 C/m 2 as a 
result of a 1.0% reduction of the in-plane lattice constant 
from the theoretical ground-state value of 7.448 bohr to 
7.374 bohr. For thicker, partially-relaxed films, the strain 
should be taken to correspond to the in-plane lattice con- 
stant measured for the film2£; this assumes that all the 
misfit dislocations responsible for the relaxation are lo- 
cated at the interface. These predictions are based on the 
assumption that the epitaxial strain strongly dominates 
other factors in determining the state of the film. 

Similarly, for strained-layer superlattices, the struc- 
ture and polarization of a perovskite-oxide layer in a su- 
perlattice will in general be significantly different from 
the corresponding bulk. The states of the component 
layers will largely be determined by the in-plane strain 
imposed by lattice matching to the other components 
and/or to the substrate, and can be obtained by refer- 
ring to the theoretical phase diagrams at the relevant 
value of misfit strain. In addition, the normal compo- 
nent of the polarization in each layer, at the common in- 
plane strain, is an important consideration in determin- 
ing the structure and polarization of the overall superlat- 
tice. If the normal polarization is discontinuous, electro- 
static energy considerations will tend to polarize the low- 
polarization layer and depolarize the high-polarization 
layer to make the normal polarization uniform through 
the sublattice^ with accompanying changes in the nor- 
mal strain of the layer. As an example, this analysis 
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has proved useful in the interpretation of first-principles 
results and experiments on the structure of partially- 
relaxed SrTi0 3 /BaTi0 3 superlattices^iSl 

First-principles information about the effects of epi- 
taxial strain on the structure and polarization of com- 
ponent layers also allows us to determine the relative 
stability of a given superlattice. It is well known that 
lattice matching of the components minimizes the elastic 
energy and thus increases the stability. In the case where 
one or more components have a nonzero normal polariza- 
tion, the electrostatic energy of the superlattice can be 
minimized by "polarization matching," where the compo- 
nents are selected to have the same normal polarization 
at the common in-plane strain. For example, according 
to our calculations, BaTi0 3 and NaNb0 3 have the same 
normal polarization of 0.34 C/m 2 at a common in-plane 
theoretical lattice constant of 7.336 bohr, corresponding 
to small compressive misfit strains of less that 0.2%, and 
thus this is a favorable combination with respect both to 
elastic and electrostatic energy. In general, however, it 
is not possible to minimize both elastic and electrostatic 
energy in this way, and a trade-off between the two is 
necessary to form the superlattice. 

It should be kept in mind, however, that for both films 
and superlattices, the assumptions that the system is in a 
single domain and that the epitaxial strain strongly dom- 
inates other factors will not be valid in all cases. Phase 
diagrams including multiple-domain states have, for ex- 
ample, been discussed in Refs. 29130,31.. Other influences 
that may be important include surface relaxation and 
reconstruction, atomic and electronic rearrangements at 
the interface, imperfectly compensated macroscopic elec- 



tric fields, deviations from stoichiometry, and the pres- 
ence of defects. 



V. SUMMARY 

We have applied the first-principles total-energy pa- 
rameterization of King-Smith and Vanderbilt^— to study 
the effects of epitaxial strain and external stress on the 
structure and properties of perovskites. We report phase 
diagrams and polarizations for the same set of eight com- 
pounds as in Ref. BaTi0 3 , SrTi0 3 , CaTi0 3 , KNb0 3 , 
NaNb0 3 , PbTi0 3 , PbZr0 3 , and BaZr0 3 . An updated 
set of parameters, computed with a comparable first- 
principles method at higher precision, are provided in 
Table IIIII The simple form of the parameterization is 
seen to be useful in reducing the computational effort 
for generating these phase diagrams relative to full first- 
principles calculations, and many features of the phase 
diagrams can be extracted and interpreted analytically. 

We have discussed the use of these results in predict- 
ing the structure and polarization of epitaxial perovskite 
films and strained layer superlattices. Additional prop- 
erties of interest, such as dielectric or piezoelectric con- 
stants can also be computed within this framework. 
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